Load Libraries

library(qiime2R) 
library(ggplot2)
library(vegan) 
## Warning: package 'vegan' was built under R version 4.1.3
## Loading required package: permute
## Warning: package 'permute' was built under R version 4.1.3
## Loading required package: lattice
## Warning: package 'lattice' was built under R version 4.1.3
## This is vegan 2.6-4
library(plyr) 
## Warning: package 'plyr' was built under R version 4.1.3
library(dplyr) 
## Warning: package 'dplyr' was built under R version 4.1.3
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:plyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(scales)
## Warning: package 'scales' was built under R version 4.1.3
library(grid) 
library(reshape2)
## Warning: package 'reshape2' was built under R version 4.1.3
library(phyloseq) 
library(picante) 
## Warning: package 'picante' was built under R version 4.1.3
## Loading required package: ape
## Warning: package 'ape' was built under R version 4.1.3
## 
## Attaching package: 'ape'
## The following object is masked from 'package:dplyr':
## 
##     where
## Loading required package: nlme
## Warning: package 'nlme' was built under R version 4.1.3
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
## 
##     collapse
library(tidyr) 
## Warning: package 'tidyr' was built under R version 4.1.3
## 
## Attaching package: 'tidyr'
## The following object is masked from 'package:reshape2':
## 
##     smiths
library(viridis) 
## Loading required package: viridisLite
## 
## Attaching package: 'viridis'
## The following object is masked from 'package:scales':
## 
##     viridis_pal
library(qiime2R) 
library(DESeq2) 
## Loading required package: S4Vectors
## Warning: package 'S4Vectors' was built under R version 4.1.3
## Loading required package: stats4
## Loading required package: BiocGenerics
## 
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:dplyr':
## 
##     combine, intersect, setdiff, union
## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
## 
##     anyDuplicated, append, as.data.frame, basename, cbind, colnames,
##     dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
##     grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
##     order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
##     rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
##     union, unique, unsplit, which.max, which.min
## 
## Attaching package: 'S4Vectors'
## The following object is masked from 'package:tidyr':
## 
##     expand
## The following objects are masked from 'package:dplyr':
## 
##     first, rename
## The following object is masked from 'package:plyr':
## 
##     rename
## The following objects are masked from 'package:base':
## 
##     expand.grid, I, unname
## Loading required package: IRanges
## 
## Attaching package: 'IRanges'
## The following object is masked from 'package:nlme':
## 
##     collapse
## The following object is masked from 'package:phyloseq':
## 
##     distance
## The following objects are masked from 'package:dplyr':
## 
##     collapse, desc, slice
## The following object is masked from 'package:plyr':
## 
##     desc
## The following object is masked from 'package:grDevices':
## 
##     windows
## Loading required package: GenomicRanges
## Loading required package: GenomeInfoDb
## Loading required package: SummarizedExperiment
## Loading required package: MatrixGenerics
## Loading required package: matrixStats
## Warning: package 'matrixStats' was built under R version 4.1.3
## 
## Attaching package: 'matrixStats'
## The following object is masked from 'package:dplyr':
## 
##     count
## The following object is masked from 'package:plyr':
## 
##     count
## 
## Attaching package: 'MatrixGenerics'
## The following objects are masked from 'package:matrixStats':
## 
##     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
##     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
##     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
##     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
##     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
##     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
##     colWeightedMeans, colWeightedMedians, colWeightedSds,
##     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
##     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
##     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
##     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
##     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
##     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
##     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
##     rowWeightedSds, rowWeightedVars
## Loading required package: Biobase
## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.
## 
## Attaching package: 'Biobase'
## The following object is masked from 'package:MatrixGenerics':
## 
##     rowMedians
## The following objects are masked from 'package:matrixStats':
## 
##     anyMissing, rowMedians
## The following object is masked from 'package:phyloseq':
## 
##     sampleNames
library(patchwork) 
library(RColorBrewer)
## Warning: package 'RColorBrewer' was built under R version 4.1.3
library(microViz) 
## microViz version 0.10.10 - Copyright (C) 2023 David Barnett
## ! Website: https://david-barnett.github.io/microViz
## v Useful?  For citation details, run: `citation("microViz")`
## x Silence? `suppressPackageStartupMessages(library(microViz))`
library(speedyseq) 
## 
## Attaching package: 'speedyseq'
## The following objects are masked from 'package:phyloseq':
## 
##     filter_taxa, plot_bar, plot_heatmap, plot_tree, psmelt, tax_glom,
##     tip_glom, transform_sample_counts
library(ComplexHeatmap) 
## ========================================
## ComplexHeatmap version 2.10.0
## Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
## Github page: https://github.com/jokergoo/ComplexHeatmap
## Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
## 
## If you use it in published research, please cite:
## Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional 
##   genomic data. Bioinformatics 2016.
## 
## The new InteractiveComplexHeatmap package can directly export static 
## complex heatmaps into an interactive Shiny app with zero effort. Have a try!
## 
## This message can be suppressed by:
##   suppressPackageStartupMessages(library(ComplexHeatmap))
## ========================================
library(ggVennDiagram) 
## Warning: package 'ggVennDiagram' was built under R version 4.1.3
library(SuperExactTest) 
## 
## Attaching package: 'SuperExactTest'
## The following objects are masked from 'package:GenomicRanges':
## 
##     intersect, union
## The following object is masked from 'package:GenomeInfoDb':
## 
##     intersect
## The following objects are masked from 'package:IRanges':
## 
##     intersect, union
## The following objects are masked from 'package:S4Vectors':
## 
##     intersect, union
## The following objects are masked from 'package:BiocGenerics':
## 
##     intersect, union
## The following objects are masked from 'package:dplyr':
## 
##     intersect, union
## The following objects are masked from 'package:base':
## 
##     intersect, union
library(nVennR)  
library(vegan) 
library(plyr) 
library(dplyr)
library(scales) 
library(grid) 
library(reshape2)
library(phyloseq) 
library(picante) 
library(tidyr) 
library(viridis) 
library(qiime2R)
library(DESeq2) 
library(patchwork) 
library(RColorBrewer)
library(speedyseq) 
library(RColorBrewer)
library(microViz) 
library(ComplexHeatmap) 
library(plotly) 
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ComplexHeatmap':
## 
##     add_heatmap
## The following object is masked from 'package:microViz':
## 
##     add_paths
## The following object is masked from 'package:IRanges':
## 
##     slice
## The following object is masked from 'package:S4Vectors':
## 
##     rename
## The following objects are masked from 'package:plyr':
## 
##     arrange, mutate, rename, summarise
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(seecolor)
library(htmlwidgets)
## Warning: package 'htmlwidgets' was built under R version 4.1.3
library(htmltools)
## Warning: package 'htmltools' was built under R version 4.1.3
set.seed(17384)

Loading Data

PRphyseq <- qza_to_phyloseq(features="S003-S015-table.qza", tree="S003-S015-rooted-tree.qza",taxonomy="S003-S015-taxonomy.qza", metadata = "S003-S015-Field-Measurements-2.txt")
rank_names(PRphyseq)
## [1] "Kingdom" "Phylum"  "Class"   "Order"   "Family"  "Genus"   "Species"
sample_variables(PRphyseq)
##  [1] "Project"                    "Sub_Location"              
##  [3] "Location"                   "Season"                    
##  [5] "Month"                      "Year"                      
##  [7] "Sample_Type"                "Sample_ID"                 
##  [9] "Day"                        "Station_ID"                
## [11] "Station_Name"               "Collected_by"              
## [13] "Time"                       "Samples_Collected_by"      
## [15] "Temperature_C"              "mm_Hg"                     
## [17] "DO_mg_L"                    "DO_LDO"                    
## [19] "Specific_Conductance_uS_cm" "Salinity_ppt"              
## [21] "pH"                         "pH_mv"                     
## [23] "IVCH_ug_L"                  "Turbidity_.Fluorometer"    
## [25] "Comments"

Removing chloroplast sequences and any contaminant sequences

## Removing chloroplast sequences and any contaminant sequences
PRphyseq <- subset_taxa(PRphyseq, Kingdom != "d__Eukaryota") 
PRphyseq <- subset_taxa(PRphyseq, Kingdom != "d__Archaea") 
PRphyseq <- subset_taxa(PRphyseq, Order != "o__Chloroplast") 
PRphyseq <- subset_taxa(PRphyseq, Order != "Chloroplast") 
PRphyseq <- subset_taxa(PRphyseq, Family != "f__Mitochondria") 
PRphyseq <- subset_taxa(PRphyseq, Family != "Mitochondria") 
PRphyseq <- subset_taxa(PRphyseq, Family != "Chloroplast")  
##Check that contaminant sequences are removed (easiest to save as data frame and search) taxtabl <- as.data.frame(tax_table(S003physeq))  #At this point, blanks and positive controls should also be removed #S003physeq = subset_samples(S003physeq, Sample != "mock")
View((tax_table(PRphyseq)))

##This will allow you to open your data frame. You can search in the top right corner of the dataframe for the sequences you wanted to remove in above code to ensure they were removed. 

Creating a subset of samples

PRphyseq <- PRphyseq %>% subset_samples(Project %in% c("S003","S015"))  ##This will create a subset (based on project column) of the samples we would like to use in the analysis. If you would like to keep all samples in the metadata file, you can ignore this step. 


PR_Normal <- PRphyseq %>% subset_samples(Sample_Type %in% c("N"))

Examine the number of reads

 readsumsdf = data.frame(nreads = sort(taxa_sums(PR_Normal), TRUE),                                           sorted = 1:ntaxa(PR_Normal), type = "ASVs") 
 readsumsdf = rbind(readsumsdf, data.frame(nreads = sort(sample_sums(PR_Normal),                          TRUE), sorted = 1:nsamples(PR_Normal), type = "Samples"))
 title = "Total number of reads" 
 p = ggplot(readsumsdf, aes(x = sorted, y = nreads)) + geom_bar(stat = "identity") 
 p + ggtitle(title) + scale_y_log10() + facet_wrap(~type, 1, scales = "free")
## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Removed 1908 rows containing missing values (`geom_bar()`).

Plot a rarefaction curve using vegan:

taxa_are_rows(PR_Normal)
## [1] TRUE
mat <- t(otu_table(PR_Normal))
class(mat) <- "matrix"
## Warning in class(mat) <- "matrix": Setting class(x) to "matrix" sets attribute
## to NULL; result will no longer be an S4 object
class(mat)
## [1] "matrix" "array"
mat <- as(t(otu_table(PR_Normal)), "matrix")
class(mat)
## [1] "matrix" "array"
raremax <- min(rowSums(mat))

system.time(rarecurve(mat, step = 100, sample = raremax, col = "blue", label = FALSE))

##    user  system elapsed 
##    4.28    0.17    6.32

Creating table for number of reads per sample

#Sum reads and sort by samples with least to greatest number of reads
samples_reads <- sort(sample_sums(PR_Normal))
#Create table
knitr::kable(samples_reads) %>% 
kableExtra::kable_styling("striped", latex_options="scale_down") %>% 
kableExtra::scroll_box(width = "100%")
x
PR109C 14
pr242a 679
PR325 699
PR164C 700
PR151C 855
PR566A 1040
PR321 1088
PR147C 1128
PR564A 1219
lsj3-1 1251
PR323 1323
PR559A 1335
PR145C 1351
PR422 1363
PR425 1372
PR162C 1450
PR141C 1461
PR421 1479
PR317 1481
PR399 1579
PR561A 1624
PR311 1634
PR420 1700
PR139C 1730
PR333 1785
PR558A 1793
PR554A 1972
PR419 2037
PR417 2075
PR160C 2191
PR107C 2201
PR556A 2311
PR334 2388
PR557A 2532
PR392 2952
PR395 3001
PR149C 3054
PR394 3055
PR332 3152
PR393 3167
PR555A 3238
PR416 3539
PR562A 3614
PR415 3728
PR328 3792
PR400 3815
PR316 4073
PR565A 4122
PR331 4337
PR402 4429
PR560A 4453
PR401 4461
PR563A 4999
PR315 5085
239 5356
18 5464
pr248a 5585
PR318 5696
PR102C 6357
PR322 6779
pr250a 7781
pr262a 8718
pr271a 9101
101 9859
pr267a 10277
220 10514
pr247a 10711
44 10755
104 10896
cs1-1 11179
cmpw-3 11765
pr295a 12285
lsj3-2 12316
cmpe-3 12325
pr249a 12448
cmp2-2 12550
142 12638
198 12984
lp1-2 13656
42 13896
cmp2-1 13943
cmpe-2 14167
cmp2-3 14339
lp1-3 14775
5 14776
218 15212
181 15472
cmpe-1 15558
cmpw-1 15630
183 16414
63 16771
lsj3-3 16802
cs1-2 19707
95 21854
lp1-1 23460
106 25062
90 25753
cs1-3 29344
cmpw-2 51074
pr293a 73188

Remove samples with <1,000 reads

#Remove samples with less than 1,000 reads
PR_Normal = subset_samples(PR_Normal, 
                            Sample_ID != "PR105C")
PR_Normal = subset_samples(PR_Normal, 
                            Sample_ID != "PR109C")
PR_Normal = subset_samples(PR_Normal, 
                           Sample_ID != "pr242a")
PR_Normal = subset_samples(PR_Normal, 
                            Sample_ID != "PR325")
PR_Normal = subset_samples(PR_Normal, 
                            Sample_ID != "PR164C")
PR_Normal = subset_samples(PR_Normal, 
                            Sample_ID != "PR151C")

Examining abundance of individual ASVs

#For each ASV, find the number of samples in which each ASV is 0, then divide by the total number of samples 
test_filter_method <- as.data.frame(rowSums(otu_table(PR_Normal) == 0)/ncol(otu_table(PR_Normal)))    
#Change the name of the column in the test_filter_method4 dataframe; this column contains the proportion of samples in which each ASV is 0  
colnames(test_filter_method) <- c('samplew0')   
#Create a bar plot 
ggplot(data = test_filter_method) + geom_bar(mapping = aes(x= samplew0)) +labs(x = "Proportion of samples in which ASV = 0", y = "# of ASV's") +theme(text = element_text(size = 18), axis.title = element_text(size = 15),panel.spacing = unit(1, "lines"),panel.border = element_rect(colour = "black", fill = NA, size = 0.5), panel.background = element_blank()) 
## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## i Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Removing singletons and doubletons

##How many singletons and doubletons are present 
#Create a dataframe that includes only the total number of times that each ASV occurs
readsumsdf_ASVs <- readsumsdf %>% 
    filter(type == "ASVs")    
#How many singletons are present?
length(which(readsumsdf_ASVs$nreads <= 0))
## [1] 1908
length(which(readsumsdf_ASVs$nreads == 1))
## [1] 7
#How many doubletons are present?
length(which(readsumsdf_ASVs$nreads == 2))
## [1] 158
#Filter out low abundance ASVs

PR_Normal = prune_taxa(taxa_sums(PR_Normal) > 2, PR_Normal)

Examine the abundance of individual ASVs again:

#For each ASV, find the number of samples in which each ASV is 0, then divide by the total number of samples 
test_filter_method_filtered <- as.data.frame(rowSums(otu_table(PR_Normal) == 0)/ncol(otu_table(PR_Normal)))   

#Change the name of the column in the test_filter_method4 dataframe; this column contains the proportion of samples in which each ASV is 0 
colnames(test_filter_method_filtered) <- c('samplew0')  

#Create a bar plot 
ggplot(data = test_filter_method_filtered) + geom_bar(mapping = aes(x= samplew0)) + 
  labs(x = "Proportion of samples in which ASV = 0", y = "# of ASV's") +   
  theme(text = element_text(size = 18), axis.title = element_text(size = 15), 
        panel.spacing = unit(1, "lines"), 
        panel.border = element_rect(colour = "black", fill = NA, size = 0.5), 
        panel.background = element_blank()) 

Label variables

sample_data(PR_Normal)$Sub_Location <- factor(sample_data(PR_Normal)$Sub_Location,              levels = c("BSJ1","BSJ1-ENR","BSJ2","BSJ2-ENR","BSJ3","BSJ3-ENR","BSJ4-ENR","BSJ5-ENR","BSJ6-ENR","BSJ7-ENR","CMP1","CMP1-ENR","CMP2","CMP2-ENR","CMPE","CMPW","CS1","CS1-ENR","CS2","CS2-ENR","CSA1-ENR","LLC1","LLC1-ENR","LP1","LP1-ENR","LP1-ENR","LP1L-ENR","LP2-ENR","LP3-ENR","LP4-ENR","LP5-ENR","LSJ1","LSJ1-ENR","LSJ2","LSJ2-ENR","LSJ3","LT1","LT1-ENR","LT2","LT2-ENR","LT3-ENR","QB1-ENR","QSA1-ENR","RPN1-ENR"),                         
labels = c("BSJ1","BSJ1-ENR","BSJ2","BSJ2-ENR","BSJ3","BSJ3-ENR","BSJ4-ENR","BSJ5-ENR","BSJ6-ENR","BSJ7-ENR","CMP1","CMP1-ENR","CMP2","CMP2-ENR","CMPE","CMPW","CS1","CS1-ENR","CS2","CS2-ENR","CSA1-ENR","LLC1","LLC1-ENR","LP1","LP1-ENR","LP1-ENR","LP1L-ENR","LP2-ENR","LP3-ENR","LP4-ENR","LP5-ENR","LSJ1","LSJ1-ENR","LSJ2","LSJ2-ENR","LSJ3","LT1","LT1-ENR","LT2","LT2-ENR","LT3-ENR","QB1-ENR","QSA1-ENR","RPN1-ENR"))

sample_data(PR_Normal)$Location <- factor(sample_data(PR_Normal)$Location,              levels = c(
  "BSJ","BSJ-ENR","CMP","CMP-ENR","CMPE","CMPW", "CS","CS-ENR", "CSA","CSA-ENR","LLC","LLC-ENR","LP","LP-ENR","LSJ","LSJ-ENR","LT","LT-ENR","QB","QB-ENR","QSA","QSA-ENR","RPN","RPN-ENR"),
labels = c("BSJ","BSJ-ENR","CMP","CMP-ENR","CMPE","CMPW", "CS","CS-ENR", "CSA","CSA-ENR","LLC","LLC-ENR","LP","LP-ENR","LSJ","LSJ-ENR","LT","LT-ENR","QB","QB-ENR","QSA","QSA-ENR","RPN","RPN-ENR"))
sample_data(PR_Normal)$Season <- factor(sample_data(PR_Normal)$Season,                                             levels = c("W", "D"), 
labels = c("W", "D"))

sample_data(PR_Normal)$Project <- factor(sample_data(PR_Normal)$Project,
                               levels = c("S003","S015"), labels = c("S003","S015"))

sample_data(PR_Normal)$Month <- factor(sample_data(PR_Normal)$Month,
                               levels = c("April","August","December","January","July","June","March","May","November","October","September"), labels = c("April","August","December","January","July","June","March","May","November","October","September"))
  
sample_data(PR_Normal)$Year <- factor(sample_data(PR_Normal)$Year,
                                     levels = c("2021","2022","2023"), labels = c("2021","2022","2023"))
  
sample_data(PR_Normal)$Sample_Type <- factor(sample_data(PR_Normal)$Sample_Type,
                                            levels = c("N","E"), labels = c("N","E"))
sample_data(PR_Normal)$Sample_ID <- factor(sample_data(PR_Normal)$Sample_ID,
                                            levels = c("PR419","PR417","PR422","PR421","PR420","PR416","PR415","90","44","42","63","5","18","pr293a", "pr271a","pr262a","pr285a","pr295a","pr267a","pr290a","PR425","PR325","PR323","PR321","PR322","PR332","PR333","PR331","PR315","PR318","PR328","PR317","PR316","PR334",    "PR311","PR563A","PR562A","PR560A","PR561A","PR565A","cmp2-3","cmpe-3","cmpw-3","cs1-3","PR556A","PR559A","PR564A","lp1-3","PR558A","PR557A","PR566A","lsj3-3","PR555A","PR554A","cmp2-2","cmpe-2","cmpw-2","cs1-2","lp1-2","lsj3-2","PR395","PR394","PR392",   "PR393","PR400","PR401","PR399","PR402","PR464A","PR465A","PR466A","PR467A","PR468A","PR469A","PR470A","PR473A","cmp2-1","PR472A","cmpe-1","cmpw-1","cs1-1","PR476A","PR477A","PR474A","PR478A","lp1-1","PR488A","PR489A","PR490A","PR491A","PR492A","PR493A","PR479A","PR480A","lsj3-1","PR481A","PR482A","PR483A","PR484A","PR485A","PR471A","Mock","Undeter","B53","B84","B59","B43","B36","B44","B107","B90","B126","B141","B128","B152","B134","B132","B140","B89","B100","B102","B60","B76","B56","blank","blank2","Mock-com","B130","B136","B138","B33","B25","B41","B81","B68","B57","B87","B91","B117","B143","B139","B154","B69","B52","B75","B99","B118","B124","B26","B34","B42","pr249a","pr248a","pr250a","pr242a","pr247a","183","218","181","198","239","220","PR160C","PR162C","PR164C","PR102C","101","PR109C","108","PR107C","106","PR143C","142","PR145C","PR151C","PR96C","95","PR147C","PR149C","PR105C","104","PR139C","PR141C"
),
labels = c("PR419","PR417","PR422","PR421","PR420","PR416","PR415","90","44","42","63","5","18","pr293a", "pr271a","pr262a","pr285a","pr295a","pr267a","pr290a","PR425","PR325","PR323","PR321","PR322","PR332","PR333","PR331","PR315","PR318","PR328","PR317","PR316","PR334",    "PR311","PR563A","PR562A","PR560A","PR561A","PR565A","cmp2-3","cmpe-3","cmpw-3","cs1-3","PR556A","PR559A","PR564A","lp1-3","PR558A","PR557A","PR566A","lsj3-3","PR555A","PR554A","cmp2-2","cmpe-2","cmpw-2","cs1-2","lp1-2","lsj3-2","PR395","PR394","PR392",   "PR393","PR400","PR401","PR399","PR402","PR464A","PR465A","PR466A","PR467A","PR468A","PR469A","PR470A","PR473A","cmp2-1","PR472A","cmpe-1","cmpw-1","cs1-1","PR476A","PR477A","PR474A","PR478A","lp1-1","PR488A","PR489A","PR490A","PR491A","PR492A","PR493A","PR479A","PR480A","lsj3-1","PR481A","PR482A","PR483A","PR484A","PR485A","PR471A","Mock","Undeter","B53","B84","B59","B43","B36","B44","B107","B90","B126","B141","B128","B152","B134","B132","B140","B89","B100","B102","B60","B76","B56","blank","blank2","Mock-com","B130","B136","B138","B33","B25","B41","B81","B68","B57","B87","B91","B117","B143","B139","B154","B69","B52","B75","B99","B118","B124","B26","B34","B42","pr249a","pr248a","pr250a","pr242a","pr247a","183","218","181","198","239","220","PR160C","PR162C","PR164C","PR102C","101","PR109C","108","PR107C","106","PR143C","142","PR145C","PR151C","PR96C","95","PR147C","PR149C","PR105C","104","PR139C","PR141C"))

Tax Fix

#View the tax_table
knitr::kable(head(tax_table(PR_Normal))) %>%   
  kableExtra::kable_styling("striped") %>%    
  kableExtra::scroll_box(width = "100%")  
Kingdom Phylum Class Order Family Genus Species
1779c60c03d948a73687f440be973f68 d__Bacteria Bacteroidota Bacteroidia Cytophagales Flammeovirgaceae Algivirga Algivirga_pacifica
25feebceed837c1dfa810cd8c0f26b41 d__Bacteria Bacteroidota Bacteroidia Cytophagales Cyclobacteriaceae NA NA
0872b535e6e15b9cd031dae6a67d6257 d__Bacteria Bacteroidota Bacteroidia Cytophagales Amoebophilaceae uncultured uncultured_bacterium
6590e42d5aeb18c338a836bc26028592 d__Bacteria Bacteroidota Bacteroidia Chitinophagales uncultured uncultured uncultured_Bacteroidetes
0f22f07d74171ba4bbe6b3977217cbdc d__Bacteria Bacteroidota Bacteroidia Bacteroidales Bacteroidetes_BD2-2 Bacteroidetes_BD2-2 uncultured_organism
d817bb1b166e13e13db6db8d56b43983 d__Bacteria Bacteroidota Bacteroidia Bacteroidales Bacteroidetes_BD2-2 Bacteroidetes_BD2-2 NA
#Now fix the labels
PR_Normal <- tax_fix(PR_Normal)

#View the relabeled tax_table
knitr::kable(head(tax_table(PR_Normal))) %>%   
  kableExtra::kable_styling("striped") %>%    
  kableExtra::scroll_box(width = "100%")
Kingdom Phylum Class Order Family Genus Species
1779c60c03d948a73687f440be973f68 d__Bacteria Bacteroidota Bacteroidia Cytophagales Flammeovirgaceae Algivirga Algivirga_pacifica
25feebceed837c1dfa810cd8c0f26b41 d__Bacteria Bacteroidota Bacteroidia Cytophagales Cyclobacteriaceae Cyclobacteriaceae Family Cyclobacteriaceae Family
0872b535e6e15b9cd031dae6a67d6257 d__Bacteria Bacteroidota Bacteroidia Cytophagales Amoebophilaceae uncultured uncultured_bacterium
6590e42d5aeb18c338a836bc26028592 d__Bacteria Bacteroidota Bacteroidia Chitinophagales uncultured uncultured uncultured_Bacteroidetes
0f22f07d74171ba4bbe6b3977217cbdc d__Bacteria Bacteroidota Bacteroidia Bacteroidales Bacteroidetes_BD2-2 Bacteroidetes_BD2-2 uncultured_organism
d817bb1b166e13e13db6db8d56b43983 d__Bacteria Bacteroidota Bacteroidia Bacteroidales Bacteroidetes_BD2-2 Bacteroidetes_BD2-2 Bacteroidetes_BD2-2 Genus

Bacterial community composition bar plots

library(seecolor) 
bcc_hex <- print_color(c("black", "#440154FF", "#450659FF","#460B5EFF","#472D7AFF","#3B518BFF", "#3A548CFF", "#38598CFF", "#365C8DFF", "#34608DFF", "#33638DFF", "#31678EFF", "#306A8EFF","#2E6E8EFF", "#2D718EFF", "#2B748EFF","#2A778EFF", "#297B8EFF","#1F958BFF", "#1F988BFF", "#1E9C89FF", "#1F9F88FF","#1FA287FF", "#21A585FF", "#23A983FF","#25AC82FF", "#29AF7FFF", "#2DB27DFF","#32B67AFF", "#37B878FF", "#3CBC74FF","#57C766FF", "#5EC962FF","#A2DA37FF","#DAE319FF", "#E4E419FF","#ECE51BFF", "#F5E61FFF","darkorchid1", "darkorchid2", "darkorchid3","#A21C9AFF", "#A62098FF", "#AB2394FF", "#AE2892FF", "#B22B8FFF", "#B6308BFF", "#BA3388FF", "#BE3885FF", "#C13B82FF", "#C53F7EFF","#C8437BFF", "#CC4678FF", "#CE4B75FF", "#D14E72FF", "#D5536FFF", "#D7566CFF", "#DA5B69FF", "#DD5E66FF", "#E06363FF","#FF6699","#F68D45FF", "#FCA537FF","#F6E726FF", "#F4ED27FF", "#00489C", "#CCCCCC", "#999999", "#A1C299","#300018"), type = "r")
## 
##  ------ c black #440154FF #450659FF #460B5EFF #472D7AFF #3B518BFF #3A548CFF #38598CFF #365C8DFF #34608DFF #33638DFF #31678EFF #306A8EFF #2E6E8EFF #2D718EFF #2B748EFF #2A778EFF #297B8EFF #1F958BFF #1F988BFF #1E9C89FF #1F9F88FF #1FA287FF #21A585FF #23A983FF #25AC82FF #29AF7FFF #2DB27DFF #32B67AFF #37B878FF #3CBC74FF #57C766FF #5EC962FF #A2DA37FF #DAE319FF #E4E419FF #ECE51BFF #F5E61FFF darkorchid1 darkorchid2 darkorchid3 #A21C9AFF #A62098FF #AB2394FF #AE2892FF #B22B8FFF #B6308BFF #BA3388FF #BE3885FF #C13B82FF #C53F7EFF #C8437BFF #CC4678FF #CE4B75FF #D14E72FF #D5536FFF #D7566CFF #DA5B69FF #DD5E66FF #E06363FF #FF6699 #F68D45FF #FCA537FF #F6E726FF #F4ED27FF #00489C #CCCCCC #999999 #A1C299 #300018 ------
## black               
## #440154FF           
## #450659FF           
## #460B5EFF           
## #472D7AFF           
## #3B518BFF           
## #3A548CFF           
## #38598CFF           
## #365C8DFF           
## #34608DFF           
## #33638DFF           
## #31678EFF           
## #306A8EFF           
## #2E6E8EFF           
## #2D718EFF           
## #2B748EFF           
## #2A778EFF           
## #297B8EFF           
## #1F958BFF           
## #1F988BFF           
## #1E9C89FF           
## #1F9F88FF           
## #1FA287FF           
## #21A585FF           
## #23A983FF           
## #25AC82FF           
## #29AF7FFF           
## #2DB27DFF           
## #32B67AFF           
## #37B878FF           
## #3CBC74FF           
## #57C766FF           
## #5EC962FF           
## #A2DA37FF           
## #DAE319FF           
## #E4E419FF           
## #ECE51BFF           
## #F5E61FFF           
## darkorchid1         
## darkorchid2         
## darkorchid3         
## #A21C9AFF           
## #A62098FF           
## #AB2394FF           
## #AE2892FF           
## #B22B8FFF           
## #B6308BFF           
## #BA3388FF           
## #BE3885FF           
## #C13B82FF           
## #C53F7EFF           
## #C8437BFF           
## #CC4678FF           
## #CE4B75FF           
## #D14E72FF           
## #D5536FFF           
## #D7566CFF           
## #DA5B69FF           
## #DD5E66FF           
## #E06363FF           
## #FF6699             
## #F68D45FF           
## #FCA537FF           
## #F6E726FF           
## #F4ED27FF           
## #00489C             
## #CCCCCC             
## #999999             
## #A1C299             
## #300018

Creating Family Barplot

##Subsetting non-enriched samples to view composition in bar plot

PR_N <- PR_Normal %>% subset_samples(Sample_Type %in% c("N"))
PRSeqR_family <- PR_N %>%  
  tax_glom(taxrank = "Family") %>%   
  transform_sample_counts(function(x) {x/sum(x)}) %>%   
  psmelt() %>%   
  group_by(Location,             
           Kingdom, Phylum, Class, Order, Family) %>%  
  filter(Abundance > 0.01) %>%    
  arrange(Class)
PRSeqR_family$Class_family <- paste(PRSeqR_family$Class, PRSeqR_family$Family, sep="_") 
PRSeqR_family_bar <- ggplot(PRSeqR_family, aes(x = Location, y = Abundance, fill = Class_family)) +   
  geom_bar(stat="identity", colour = "black", size=0.3, position = "fill") + scale_fill_manual(values = c("black",                                                                 "#440154FF", "#450659FF","#460B5EFF",                                                                "#472D7AFF",                                                                 "#3B518BFF", "#3A548CFF", "#38598CFF", "#365C8DFF",                                 "#34608DFF", "#33638DFF", "#31678EFF", "#306A8EFF",                                "#2E6E8EFF", "#2D718EFF", "#2B748EFF",                                "#2A778EFF", "#297B8EFF",                                                                "#1F958BFF", "#1F988BFF", "#1E9C89FF", "#1F9F88FF",                                "#1FA287FF", "#21A585FF", "#23A983FF",                                "#25AC82FF", "#29AF7FFF", "#2DB27DFF",                                "#32B67AFF", "#37B878FF", "#3CBC74FF",                                                                "#57C766FF", "#5EC962FF",                                                                "#A2DA37FF",                                                                "#DAE319FF", "#E4E419FF",  "#ECE51BFF", "#F5E61FFF",                                                                "darkorchid1", "darkorchid2", "darkorchid3",                                                                 "#A21C9AFF", "#A62098FF", "#AB2394FF", "#AE2892FF", "#B22B8FFF", "#B6308BFF", "#BA3388FF", "#BE3885FF", "#C13B82FF", "#C53F7EFF",                                "#C8437BFF", "#CC4678FF", "#CE4B75FF", "#D14E72FF", "#D5536FFF", "#D7566CFF", "#DA5B69FF", "#DD5E66FF", "#E06363FF","#00CCCC", "#006600",                                                                "#FF6699",                                                                "#F68D45FF", "#FCA537FF",                                                                "#F6E726FF", "#F4ED27FF",                                                                                                 "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018",                                                                "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018",                                "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018",                                "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018","#191970","#353935","#F0FFFF","#89CFF0","#0000FF","#7393B3","#088F8F","#0096FF","#5F9EA0","#0047AB","#6495ED","#00FFFF","#00008B","#6F8FAF","#1434A4","#7DF9FF","#6082B6","#00A36C","#3F00FF","#5D3FD3","#ADD8E6","#191970","#000080","#1F51FF","#A7C7E7","#CCCCFF","#B6D0E2","#96DED1","#4169E1","#0F52BA","#9FE2BF","#87CEEB","#4682B4","#008080","#40E0D0","#0437F2","#40B5AD","#0818A8","#AAFF00","#5F9EA0","#097969","#AFE1AF","#DFFF00","#E4D00A","#00FFFF","#023020","#50C878","#5F8575","#4F7942","#228B22","#7CFC00","#F2D2BD","#FFAC1C","#CD7F32","#DAA06D","#CC5500","#E97451","#E3963E","#F28C28","#D27D2D","#B87333","#FF7F50","#F88379","#DA70D6","#C3B1E1","#CCCCFF","#673147","#A95C68","#800080","#51414F","#953553","#D8BFD8","#630330","#7F00FF","#722F37","#BDB5D5","#FCF55F","#FFFFF0","#F8DE7E","#F0E68C","#FAFA33","#FBEC5D","#F4BB44","#FFDB58")) +  
  guides(fill = guide_legend(reverse = FALSE, keywidth = 1, keyheight = 1)) +  
  ylab("Relative Abundance (Family > 1%) \n") + xlab("Location") +   
  theme_minimal()+theme(panel.border = element_rect(fill=NA, colour = "black"), panel.grid.minor = element_blank(),                         
panel.grid.major = element_blank(), axis.text.x  =element_text(angle = 90, vjust = 0.5, hjust=1, size=8, colour="black"),                          
axis.text.y = element_text(size=8, colour="black"), plot.title = element_text(hjust = 0.5),     
axis.ticks.x = element_line(colour="#000000", size=0.1), axis.ticks.y = element_line(colour="#000000", size=0.1),                        
legend.text = element_text(size = 7))
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## i Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## i Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
print(PRSeqR_family_bar)

p <- ggplotly(PRSeqR_family_bar)
htmlwidgets::saveWidget(p, "PRSeqR_family_bar.html")

htmltools::tags$iframe(
  src=file.path(getwd(), "PRSeqR_family_bar.html"),
  width="100%",
  height="600",
  scrolling="no",
  seamless="seamless",
  frameBorder="0")

Creating Species Barplot

PRSeqR_Species <- PR_N %>%  
  tax_glom(taxrank = "Species") %>%   
  transform_sample_counts(function(x) {x/sum(x)}) %>%   
  psmelt() %>%   
  group_by(Location,             
           Kingdom, Phylum, Class, Order, Family,Genus,Species) %>%  
  filter(Abundance > 0.05) %>%    
  arrange(Species)
PRSeqR_Species_bar <- ggplot(PRSeqR_Species, aes(x = Location, y = Abundance, fill = Species)) +   
  geom_bar(stat="identity", colour = "black", size=0.3, position = "fill") + scale_fill_manual(values = c("black",                                                                 "#440154FF", "#450659FF","#460B5EFF",                                                                "#472D7AFF",                                                                 "#3B518BFF", "#3A548CFF", "#38598CFF", "#365C8DFF",                                 "#34608DFF", "#33638DFF", "#31678EFF", "#306A8EFF",                                "#2E6E8EFF", "#2D718EFF", "#2B748EFF",                                "#2A778EFF", "#297B8EFF",                                                                "#1F958BFF", "#1F988BFF", "#1E9C89FF", "#1F9F88FF",                                "#1FA287FF", "#21A585FF", "#23A983FF",                                "#25AC82FF", "#29AF7FFF", "#2DB27DFF",                                "#32B67AFF", "#37B878FF", "#3CBC74FF",                                                                "#57C766FF", "#5EC962FF",                                                                "#A2DA37FF",                                                                "#DAE319FF", "#E4E419FF",  "#ECE51BFF", "#F5E61FFF",                                                                "darkorchid1", "darkorchid2", "darkorchid3",                                                                 "#A21C9AFF", "#A62098FF", "#AB2394FF", "#AE2892FF", "#B22B8FFF", "#B6308BFF", "#BA3388FF", "#BE3885FF", "#C13B82FF", "#C53F7EFF",                                "#C8437BFF", "#CC4678FF", "#CE4B75FF", "#D14E72FF", "#D5536FFF", "#D7566CFF", "#DA5B69FF", "#DD5E66FF", "#E06363FF","#00CCCC", "#006600",                                                                "#FF6699",                                                                "#F68D45FF", "#FCA537FF",                                                                "#F6E726FF", "#F4ED27FF",                                                                                                 "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018",                                                                "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018",                                "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018",                                "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018","#191970","#353935","#F0FFFF","#89CFF0","#0000FF","#7393B3","#088F8F","#0096FF","#5F9EA0","#0047AB","#6495ED","#00FFFF","#00008B","#6F8FAF","#1434A4","#7DF9FF","#6082B6","#00A36C","#3F00FF","#5D3FD3","#ADD8E6","#191970","#000080","#1F51FF","#A7C7E7","#CCCCFF","#B6D0E2","#96DED1","#4169E1","#0F52BA","#9FE2BF","#87CEEB","#4682B4","#008080","#40E0D0","#0437F2","#40B5AD","#0818A8","#AAFF00","#5F9EA0","#097969","#AFE1AF","#DFFF00","#E4D00A","#00FFFF","#023020","#50C878","#5F8575","#4F7942","#228B22","#7CFC00","#F2D2BD","#FFAC1C","#CD7F32","#DAA06D","#CC5500","#E97451","#E3963E","#F28C28","#D27D2D","#B87333","#FF7F50","#F88379","#DA70D6","#C3B1E1","#CCCCFF","#673147","#A95C68","#800080","#51414F","#953553","#D8BFD8","#630330","#7F00FF","#722F37","#BDB5D5","#FCF55F","#FFFFF0","#F8DE7E","#F0E68C","#FAFA33","#FBEC5D","#F4BB44","#FFDB58")) +  
  guides(fill = guide_legend(reverse = FALSE, keywidth = 1, keyheight = 1)) +  
  ylab("Relative Abundance (Species > 5%) \n") + xlab("Sample ID") +   
  theme_minimal()+theme(panel.border = element_rect(fill=NA, colour = "black"), panel.grid.minor = element_blank(),                         
panel.grid.major = element_blank(), axis.text.x  =element_text(angle = 90, vjust = 0.5, hjust=1, size=8, colour="black"),                          
axis.text.y = element_text(size=8, colour="black"), plot.title = element_text(hjust = 0.5),     
axis.ticks.x = element_line(colour="#000000", size=0.1), axis.ticks.y = element_line(colour="#000000", size=0.1),                        
legend.text = element_text(size = 7))
p2 <- ggplotly(PRSeqR_Species_bar)
print(PRSeqR_Species_bar)

htmlwidgets::saveWidget(p, "PRSeqR_Species_bar.html")

htmltools::tags$iframe(
  src=file.path(getwd(), "PRSeqR_Species_bar.html"),
  width="100%",
  height="600",
  scrolling="no",
  seamless="seamless",
  frameBorder="0")

Creating Season vs Location Barplot with non-enriched samples

PRSeqR_family_Season <- PR_N %>%    
  tax_glom(taxrank = "Family") %>%     
  transform_sample_counts(function(x) {x/sum(x)}) %>%   
  psmelt() %>%    
  group_by(Season,             
           Kingdom, Phylum, Class, Order, Family)%>%   
  dplyr::summarize(Mean =                      
                     mean(Abundance, na.rm=TRUE)) %>%    
  filter(Mean > 0.01) %>%                                
  arrange(Class)  
## `summarise()` has grouped output by 'Season', 'Kingdom', 'Phylum', 'Class',
## 'Order'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'Season', 'Kingdom', 'Phylum', 'Class',
## 'Order'. You can override using the `.groups` argument.
## Creaitng the plot
PRSeqR_family_Season$Class_family <- paste(PRSeqR_family_Season$Class, PRSeqR_family_Season$Family, sep="_")  
PRSeqR_Season_family_bar <- ggplot(PRSeqR_family_Season, aes(x = Season, y = Mean, fill = Class_family)) +   
  geom_bar(stat="identity", colour = "black", size=0.3, position = "fill") +   scale_fill_manual(values = c("black", "#440154FF", "#450659FF","#460B5EFF","#472D7AFF", "#3B518BFF", "#3A548CFF", "#38598CFF", "#365C8DFF","#34608DFF", "#33638DFF", "#31678EFF", "#306A8EFF","#2E6E8EFF", "#2D718EFF", "#2B748EFF","#2A778EFF", "#297B8EFF","#1F958BFF", "#1F988BFF", "#1E9C89FF", "#1F9F88FF", "#1FA287FF", "#21A585FF", "#23A983FF","#25AC82FF", "#29AF7FFF", "#2DB27DFF","#32B67AFF", "#37B878FF", "#3CBC74FF","#57C766FF", "#5EC962FF","#A2DA37FF", "#DAE319FF", "#E4E419FF",  "#ECE51BFF", "#F5E61FFF","darkorchid1", "darkorchid2", "darkorchid3", "#A21C9AFF", "#A62098FF", "#AB2394FF", "#AE2892FF", "#B22B8FFF", "#B6308BFF", "#BA3388FF", "#BE3885FF", "#C13B82FF", "#C53F7EFF","#C8437BFF", "#CC4678FF", "#CE4B75FF", "#D14E72FF", "#D5536FFF", "#D7566CFF", "#DA5B69FF", "#DD5E66FF", "#E06363FF","#FF6699","#F68D45FF", "#FCA537FF","#F6E726FF", "#F4ED27FF", "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018","#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018","#00489C", "#CCCCCC", "#999999","#999999", "#A1C299", "#300018","#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018")) +      guides(fill = guide_legend(reverse = FALSE, keywidth = 1, keyheight = 1)) +   ylab("Relative Abundance (Family > 1%) \n") + xlab("Sample ID") +   theme_minimal()+theme(panel.border = element_rect(fill=NA, colour = "black"), panel.grid.minor = element_blank(),                          panel.grid.major = element_blank(), axis.text.x  = element_text(angle = 90, vjust = 0.5, hjust=1, size=8, colour="black"),                          axis.text.y = element_text(size=8, colour="black"), plot.title = element_text(hjust = 0.5),                          axis.ticks.x = element_line(colour="#000000", size=0.1), axis.ticks.y = element_line(colour="#000000", size=0.1),                         legend.text = element_text(size = 7))    
print(PRSeqR_Season_family_bar)

p3 <- ggplotly(PRSeqR_Season_family_bar)

p
htmlwidgets::saveWidget(p, "Season_family.html")

htmltools::tags$iframe(
  src=file.path(getwd(), "Season_family.html"),
  width="100%",
  height="600",
  scrolling="no",
  seamless="seamless",
  frameBorder="0")

Creating Season vs. Location barplot with only normal samples

PRSeqR_family_Location <- PR_Normal %>%   
  tax_glom(taxrank = "Family") %>%                    
  #Agglomerate at family level   
  transform_sample_counts(function(x) {x/sum(x)}) %>%
  #Transform to rel.abundance  
  psmelt() %>%                                    
  # Melt to long format  
  group_by(Season, Location, Kingdom, Phylum, Class, Order, Family) %>%   
  dplyr::summarize(Mean =  mean(Abundance, na.rm=TRUE)) %>%                              #Calculate average  
  filter(Mean > 0.01) %>%                         
  #Filter arrange
  arrange(Class) 
## `summarise()` has grouped output by 'Season', 'Location', 'Kingdom', 'Phylum',
## 'Class', 'Order'. You can override using the `.groups` argument.
#Creating a new column that combines both the Class and Family level information #so that Family level identifiers that are not unique to one Class (such as Ambiguous_taxa) aren't merged into a single large category when graphed 
PRSeqR_family_Location$Class_family <- paste(PRSeqR_family_Location$Class, PRSeqR_family_Location$Family, sep="_")                                                 #Creating plot   
PRSeqR_Location_family_bar <- ggplot(PRSeqR_family_Location, aes(x = Location, y = Mean, fill = Class_family)) +  geom_bar(stat="identity", colour = "black", size=0.3, position = "fill") +   scale_fill_manual(values = c("black",                                                                 "#440154FF", "#450659FF","#460B5EFF",                                                                "#472D7AFF",                                                                 "#3B518BFF", "#3A548CFF", "#38598CFF", "#365C8DFF",                                 "#34608DFF", "#33638DFF", "#31678EFF", "#306A8EFF",                                "#2E6E8EFF", "#2D718EFF", "#2B748EFF",                                "#2A778EFF", "#297B8EFF",                                                                "#1F958BFF", "#1F988BFF", "#1E9C89FF", "#1F9F88FF",                                "#1FA287FF", "#21A585FF", "#23A983FF",                                "#25AC82FF", "#29AF7FFF", "#2DB27DFF",                                "#32B67AFF", "#37B878FF", "#3CBC74FF",                                                                "#57C766FF", "#5EC962FF",                                                                "#A2DA37FF",                                                                "#DAE319FF", "#E4E419FF",  "#ECE51BFF", "#F5E61FFF",                                                                "darkorchid1", "darkorchid2", "darkorchid3",                                                                 "#A21C9AFF", "#A62098FF", "#AB2394FF", "#AE2892FF", "#B22B8FFF", "#B6308BFF", "#BA3388FF", "#BE3885FF", "#C13B82FF", "#C53F7EFF",                                "#C8437BFF", "#CC4678FF", "#CE4B75FF", "#D14E72FF", "#D5536FFF", "#D7566CFF", "#DA5B69FF", "#DD5E66FF", "#E06363FF","#00CCCC", "#006600",                                                                "#FF6699",                                                                "#F68D45FF", "#FCA537FF",                                                                "#F6E726FF", "#F4ED27FF",                                                                                                 "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018",                                                                "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018",                                "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018",                                "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018","#191970","#353935","#F0FFFF","#89CFF0","#0000FF","#7393B3","#088F8F","#0096FF","#5F9EA0","#0047AB","#6495ED","#00FFFF","#00008B","#6F8FAF","#1434A4","#7DF9FF","#6082B6","#00A36C","#3F00FF","#5D3FD3","#ADD8E6","#191970","#000080","#1F51FF","#A7C7E7","#CCCCFF","#B6D0E2","#96DED1","#4169E1","#0F52BA","#9FE2BF","#87CEEB","#4682B4","#008080","#40E0D0","#0437F2","#40B5AD","#0818A8","#AAFF00","#5F9EA0","#097969","#AFE1AF","#DFFF00","#E4D00A","#00FFFF","#023020","#50C878","#5F8575","#4F7942","#228B22","#7CFC00","#F2D2BD","#FFAC1C","#CD7F32","#DAA06D","#CC5500","#E97451","#E3963E","#F28C28","#D27D2D","#B87333","#FF7F50","#F88379","#DA70D6","#C3B1E1","#CCCCFF","#673147","#A95C68","#800080","#51414F","#953553","#D8BFD8","#630330","#7F00FF","#722F37","#BDB5D5","#FCF55F","#FFFFF0","#F8DE7E","#F0E68C","#FAFA33","#FBEC5D","#F4BB44","#FFDB58")) +      guides(fill = guide_legend(reverse = FALSE, keywidth = 1, keyheight = 1)) +   ylab("Relative Abundance (Family > 1%) \n") + xlab("Location") +   theme_minimal()+theme(panel.border = element_rect(fill=NA, colour = "black"),    panel.grid.minor = element_blank(),   panel.grid.major = element_blank(),    axis.text.x  = element_text(angle = 90, vjust = 0.5, hjust=1, size=8,                                colour="black"),    axis.text.y = element_text(size=8, colour="black"),    plot.title = element_text(hjust = 0.5),    axis.ticks.x = element_line(colour="#000000", size=0.1),    axis.ticks.y = element_line(colour="#000000", size=0.1),   legend.text = element_text(size = 7)) +   facet_grid(. ~ Location, margins = TRUE, scale="free")

print(PRSeqR_Location_family_bar)

p4 <- ggplotly(PRSeqR_Location_family_bar)
htmlwidgets::saveWidget(p,"PRSeqR_Location_family_bar")
htmltools::tags$iframe(
  src=file.path(getwd(),"PRSeqR_Location_family_bar"),
  width="100%",
  height="600",
  scrolling="no",
  seamless="seamless",
  frameBorder="0")

Creating location vs. year barplot

PR_Normal <- PRphyseq %>% subset_samples(Sample_Type %in% c("N"))

PR_2022 <- PRphyseq %>% subset_samples(Year %in% c("2022"))

##Creating a subset of samples for only normal samples.This will exclude the enriched samples
PRSeqR_family_Year <- PRphyseq %>%    
  tax_glom(taxrank = "Family") %>%     
  transform_sample_counts(function(x) {x/sum(x)}) %>%   
  psmelt() %>%    
  group_by(Year,Season, Kingdom, Phylum, Class, Order, Family)%>%   
  dplyr::summarize(Mean =                      
                     mean(Abundance, na.rm=TRUE)) %>%    
  filter(Mean > 0.01) %>%                                
  arrange(Class)  
## `summarise()` has grouped output by 'Year', 'Season', 'Kingdom', 'Phylum',
## 'Class', 'Order'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'Season', 'Kingdom', 'Phylum', 'Class',
## 'Order'. You can override using the `.groups` argument.
#Creating a new column that combines both the Class and Family level information #so that Family level identifiers that are not unique to one Class (such as Ambiguous_taxa) aren't merged into a single large category when graphed 
PRSeqR_family_Year$Class_family <- paste(PRSeqR_family_Year$Class, PRSeqR_family_Year$Family, sep="_")                                                 #Creating plot   
PRSeqR_Year_family_bar <- ggplot(PRSeqR_family_Year, aes(x = Year, y = Mean, fill = Class_family)) +  geom_bar(stat="identity", colour = "black", size=0.3, position = "fill") +   scale_fill_manual(values = c("black",                                                                 "#440154FF", "#450659FF","#460B5EFF",                                                                "#472D7AFF",                                                                 "#3B518BFF", "#3A548CFF", "#38598CFF", "#365C8DFF",                                 "#34608DFF", "#33638DFF", "#31678EFF", "#306A8EFF",                                "#2E6E8EFF", "#2D718EFF", "#2B748EFF",                                "#2A778EFF", "#297B8EFF",                                                                "#1F958BFF", "#1F988BFF", "#1E9C89FF", "#1F9F88FF",                                "#1FA287FF", "#21A585FF", "#23A983FF",                                "#25AC82FF", "#29AF7FFF", "#2DB27DFF",                                "#32B67AFF", "#37B878FF", "#3CBC74FF",                                                                "#57C766FF", "#5EC962FF",                                                                "#A2DA37FF",                                                                "#DAE319FF", "#E4E419FF",  "#ECE51BFF", "#F5E61FFF",                                                                "darkorchid1", "darkorchid2", "darkorchid3",                                                                 "#A21C9AFF", "#A62098FF", "#AB2394FF", "#AE2892FF", "#B22B8FFF", "#B6308BFF", "#BA3388FF", "#BE3885FF", "#C13B82FF", "#C53F7EFF",                                "#C8437BFF", "#CC4678FF", "#CE4B75FF", "#D14E72FF", "#D5536FFF", "#D7566CFF", "#DA5B69FF", "#DD5E66FF", "#E06363FF","#00CCCC", "#006600",                                                                "#FF6699",                                                                "#F68D45FF", "#FCA537FF",                                                                "#F6E726FF", "#F4ED27FF",                                                                                                 "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018",                                                                "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018",                                "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018",                                "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018","#191970","#353935","#F0FFFF","#89CFF0","#0000FF","#7393B3","#088F8F","#0096FF","#5F9EA0","#0047AB","#6495ED","#00FFFF","#00008B","#6F8FAF","#1434A4","#7DF9FF","#6082B6","#00A36C","#3F00FF","#5D3FD3","#ADD8E6","#191970","#000080","#1F51FF","#A7C7E7","#CCCCFF","#B6D0E2","#96DED1","#4169E1","#0F52BA","#9FE2BF","#87CEEB","#4682B4","#008080","#40E0D0","#0437F2","#40B5AD","#0818A8","#AAFF00","#5F9EA0","#097969","#AFE1AF","#DFFF00","#E4D00A","#00FFFF","#023020","#50C878","#5F8575","#4F7942","#228B22","#7CFC00","#F2D2BD","#FFAC1C","#CD7F32","#DAA06D","#CC5500","#E97451","#E3963E","#F28C28","#D27D2D","#B87333","#FF7F50","#F88379","#DA70D6","#C3B1E1","#CCCCFF","#673147","#A95C68","#800080","#51414F","#953553","#D8BFD8","#630330","#7F00FF","#722F37","#BDB5D5","#FCF55F","#FFFFF0","#F8DE7E","#F0E68C","#FAFA33","#FBEC5D","#F4BB44","#FFDB58")) +      guides(fill = guide_legend(reverse = FALSE, keywidth = 1, keyheight = 1)) +   ylab("Relative Abundance (Family > 1%) \n") + xlab("Year") +   theme_minimal()+theme(panel.border = element_rect(fill=NA, colour = "black"),    panel.grid.minor = element_blank(),   panel.grid.major = element_blank(),    axis.text.x  = element_text(angle = 90, vjust = 0.5, hjust=1, size=8,                                colour="black"),    axis.text.y = element_text(size=8, colour="black"),    plot.title = element_text(hjust = 0.5),    axis.ticks.x = element_line(colour="#000000", size=0.1),    axis.ticks.y = element_line(colour="#000000", size=0.1),   legend.text = element_text(size = 7)) +   facet_grid(. ~ Season, margins = TRUE, scale="free")

print(PRSeqR_Year_family_bar)

p5 <- ggplotly(PRSeqR_Year_family_bar)

Note that 2023 is only enriched samples so this year cannot be accurately depicted and 2022 is missing from the wet season.

*Now that we know there is an error in the 2022 data when attempting to generate the 2022 NMDS plot, we will be creating a family plot for the 2022 samples to determine what is too similar - this is what we think is preventing the NMDS plot from being generated.

PRSeqR_2022_Sample <- PR_2022 %>%
  tax_glom(taxrank = "Family") %>%
  transform_sample_counts(function(x) {x/sum(x)}) %>%
  psmelt() %>%
  group_by(Location,Season,Kingdom,Phylum,Class,Order,Family,Genus,Species)%>%   
  dplyr::summarize(Mean =          
 mean(Abundance, na.rm=TRUE)) %>%    
  filter(Mean > 0.00) %>%                                
  arrange(Class)  
## `summarise()` has grouped output by 'Location', 'Season', 'Kingdom', 'Phylum',
## 'Class', 'Order', 'Family', 'Genus'. You can override using the `.groups`
## argument.
##`summarise()` has grouped output by 'Season', 'Kingdom', 'Phylum', 'Class',
## 'Order'. You can override using the `.groups` argument
#Creating a new column that combines both the Class and Family level information #so that Family level identifiers that are not unique to one Class (such as Ambiguous_taxa) aren't merged into a single large category when graphed 
PRSeqR_2022_Sample$Family <- paste(PRSeqR_2022_Sample$Class, PRSeqR_2022_Sample$Family, sep="_")                                                 #Creating plot   
PRSeqR_2022_Sample_bar <- ggplot(PRSeqR_2022_Sample, aes(x = Location, y = Mean, fill = Family)) +  geom_bar(stat="identity", colour = "black", size=0.3, position = "fill") +   scale_fill_manual(values = c("black",                                                                 "#440154FF", "#450659FF","#460B5EFF",                                                                "#472D7AFF",                                                                 "#3B518BFF", "#3A548CFF", "#38598CFF", "#365C8DFF",                                 "#34608DFF", "#33638DFF", "#31678EFF", "#306A8EFF",                                "#2E6E8EFF", "#2D718EFF", "#2B748EFF",                                "#2A778EFF", "#297B8EFF",                                                                "#1F958BFF", "#1F988BFF", "#1E9C89FF", "#1F9F88FF",                                "#1FA287FF", "#21A585FF", "#23A983FF",                                "#25AC82FF", "#29AF7FFF", "#2DB27DFF",                                "#32B67AFF", "#37B878FF", "#3CBC74FF",                                                                "#57C766FF", "#5EC962FF",                                                                "#A2DA37FF",                                                                "#DAE319FF", "#E4E419FF",  "#ECE51BFF", "#F5E61FFF",                                                                "darkorchid1", "darkorchid2", "darkorchid3",                                                                 "#A21C9AFF", "#A62098FF", "#AB2394FF", "#AE2892FF", "#B22B8FFF", "#B6308BFF", "#BA3388FF", "#BE3885FF", "#C13B82FF", "#C53F7EFF",                                "#C8437BFF", "#CC4678FF", "#CE4B75FF", "#D14E72FF", "#D5536FFF", "#D7566CFF", "#DA5B69FF", "#DD5E66FF", "#E06363FF","#00CCCC", "#006600",                                                                "#FF6699",                                                                "#F68D45FF", "#FCA537FF",                                                                "#F6E726FF", "#F4ED27FF",                                                                                                 "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018",                                                                "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018",                                "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018",                                "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018","#191970","#353935","#F0FFFF","#89CFF0","#0000FF","#7393B3","#088F8F","#0096FF","#5F9EA0","#0047AB","#6495ED","#00FFFF","#00008B","#6F8FAF","#1434A4","#7DF9FF","#6082B6","#00A36C","#3F00FF","#5D3FD3","#ADD8E6","#191970","#000080","#1F51FF","#A7C7E7","#CCCCFF","#B6D0E2","#96DED1","#4169E1","#0F52BA","#9FE2BF","#87CEEB","#4682B4","#008080","#40E0D0","#0437F2","#40B5AD","#0818A8","#AAFF00","#5F9EA0","#097969","#AFE1AF","#DFFF00","#E4D00A","#00FFFF","#023020","#50C878","#5F8575","#4F7942","#228B22","#7CFC00","#F2D2BD","#FFAC1C","#CD7F32","#DAA06D","#CC5500","#E97451","#E3963E","#F28C28","#D27D2D","#B87333","#FF7F50","#F88379","#DA70D6","#C3B1E1","#CCCCFF","#673147","#A95C68","#800080","#51414F","#953553","#D8BFD8","#630330","#7F00FF","#722F37","#BDB5D5","#FCF55F","#FFFFF0","#F8DE7E","#F0E68C","#FAFA33","#FBEC5D","#F4BB44","#FFDB58")) +      guides(fill = guide_legend(reverse = FALSE, keywidth = 1, keyheight = 1)) +   ylab("Relative Abundance (Family > 1%) \n") + xlab("Location") +   theme_minimal()+theme(panel.border = element_rect(fill=NA, colour = "black"),    panel.grid.minor = element_blank(),   panel.grid.major = element_blank(),    axis.text.x  = element_text(angle = 90, vjust = 0.5, hjust=1, size=8,                                colour="black"),    axis.text.y = element_text(size=8, colour="black"),    plot.title = element_text(hjust = 0.5),    axis.ticks.x = element_line(colour="#000000", size=0.1),    axis.ticks.y = element_line(colour="#000000", size=0.1),   legend.text = element_text(size = 7))

print(PRSeqR_2022_Sample_bar)

 p6 <- ggplotly(PRSeqR_2022_Sample_bar)

This is showing the samples averaged together, but there are multiple samples from each site. We are going to plot based on sample name to give us individual bars for each sample to better understand the similarities between each other.

PRSeqR_2022_Sample_ID <- PR_2022 %>%    
  tax_glom(taxrank = "Family") %>%     
  transform_sample_counts(function(x) {x/sum(x)}) %>%   
  psmelt() %>%    
  group_by(Sample_ID,Location,Season,Kingdom,Phylum,Class,Order,Family,Genus,Species)%>%   
  dplyr::summarize(Mean =                      
                     mean(Abundance, na.rm=TRUE)) %>%    
  filter(Mean > 0.00) %>%                                
  arrange(Class)  
## `summarise()` has grouped output by 'Sample_ID', 'Location', 'Season',
## 'Kingdom', 'Phylum', 'Class', 'Order', 'Family', 'Genus'. You can override
## using the `.groups` argument.
## `summarise()` has grouped output by 'Season', 'Kingdom', 'Phylum', 'Class',
## 'Order'. You can override using the `.groups` argument
#Creating a new column that combines both the Class and Family level information #so that Family level identifiers that are not unique to one Class (such as Ambiguous_taxa) aren't merged into a single large category when graphed 
PRSeqR_2022_Sample_ID$Family <- paste(PRSeqR_2022_Sample_ID$Class, PRSeqR_2022_Sample_ID$Family, sep="_")                                                 #Creating plot   
PRSeqR_2022_Sample_ID_bar <- ggplot(PRSeqR_2022_Sample_ID, aes(x = Sample_ID, y = Mean, fill = Family)) +  geom_bar(stat="identity", colour = "black", size=0.3, position = "fill") +   scale_fill_manual(values = c("black",                                                                 "#440154FF", "#450659FF","#460B5EFF",                                                                "#472D7AFF",                                                                 "#3B518BFF", "#3A548CFF", "#38598CFF", "#365C8DFF",                                 "#34608DFF", "#33638DFF", "#31678EFF", "#306A8EFF",                                "#2E6E8EFF", "#2D718EFF", "#2B748EFF",                                "#2A778EFF", "#297B8EFF",                                                                "#1F958BFF", "#1F988BFF", "#1E9C89FF", "#1F9F88FF",                                "#1FA287FF", "#21A585FF", "#23A983FF",                                "#25AC82FF", "#29AF7FFF", "#2DB27DFF",                                "#32B67AFF", "#37B878FF", "#3CBC74FF",                                                                "#57C766FF", "#5EC962FF",                                                                "#A2DA37FF",                                                                "#DAE319FF", "#E4E419FF",  "#ECE51BFF", "#F5E61FFF",                                                                "darkorchid1", "darkorchid2", "darkorchid3",                                                                 "#A21C9AFF", "#A62098FF", "#AB2394FF", "#AE2892FF", "#B22B8FFF", "#B6308BFF", "#BA3388FF", "#BE3885FF", "#C13B82FF", "#C53F7EFF",                                "#C8437BFF", "#CC4678FF", "#CE4B75FF", "#D14E72FF", "#D5536FFF", "#D7566CFF", "#DA5B69FF", "#DD5E66FF", "#E06363FF","#00CCCC", "#006600",                                                                "#FF6699",                                                                "#F68D45FF", "#FCA537FF",                                                                "#F6E726FF", "#F4ED27FF",                                                                                                 "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018",                                                                "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018",                                "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018",                                "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018","#191970","#353935","#F0FFFF","#89CFF0","#0000FF","#7393B3","#088F8F","#0096FF","#5F9EA0","#0047AB","#6495ED","#00FFFF","#00008B","#6F8FAF","#1434A4","#7DF9FF","#6082B6","#00A36C","#3F00FF","#5D3FD3","#ADD8E6","#191970","#000080","#1F51FF","#A7C7E7","#CCCCFF","#B6D0E2","#96DED1","#4169E1","#0F52BA","#9FE2BF","#87CEEB","#4682B4","#008080","#40E0D0","#0437F2","#40B5AD","#0818A8","#AAFF00","#5F9EA0","#097969","#AFE1AF","#DFFF00","#E4D00A","#00FFFF","#023020","#50C878","#5F8575","#4F7942","#228B22","#7CFC00","#F2D2BD","#FFAC1C","#CD7F32","#DAA06D","#CC5500","#E97451","#E3963E","#F28C28","#D27D2D","#B87333","#FF7F50","#F88379","#DA70D6","#C3B1E1","#CCCCFF","#673147","#A95C68","#800080","#51414F","#953553","#D8BFD8","#630330","#7F00FF","#722F37","#BDB5D5","#FCF55F","#FFFFF0","#F8DE7E","#F0E68C","#FAFA33","#FBEC5D","#F4BB44","#FFDB58")) +      guides(fill = guide_legend(reverse = FALSE, keywidth = 1, keyheight = 1)) +   ylab("Relative Abundance (Family > 1%) \n") + xlab("Sample_ID") +   theme_minimal()+theme(panel.border = element_rect(fill=NA, colour = "black"),    panel.grid.minor = element_blank(),   panel.grid.major = element_blank(),    axis.text.x  = element_text(angle = 90, vjust = 0.5, hjust=1, size=8,                                colour="black"),    axis.text.y = element_text(size=8, colour="black"),    plot.title = element_text(hjust = 0.5),    axis.ticks.x = element_line(colour="#000000", size=0.1),    axis.ticks.y = element_line(colour="#000000", size=0.1),   legend.text = element_text(size = 7))

print(PRSeqR_2022_Sample_ID_bar)

p7 <- ggplotly(PRSeqR_2022_Sample_ID_bar)
htmlwidgets::saveWidget(p,"PRSeqR_2022_Sample_ID_bar")
htmltools::tags$iframe(
  src=file.path(getwd(),"PRSeqR_2022_Sample_ID_bar"),
  width="100%",
  height="600",
  scrolling="no",
  seamless="seamless",
  frameBorder="0")

Creating barplots for enriched samples

#PR_ENR <- PRphyseq %>% subset_samples(Sample_Type %in% c("E"))

##Creating a subset of samples for only enriched samples.This will exclude the normal samples
#PRSeqR_family_Location <- PR_ENR %>%    
#  tax_glom(taxrank = "Species") %>%     
#  transform_sample_counts(function(x) {x/sum(x)}) %>%   
#  psmelt() %>%    
#  group_by(Location,Season,Kingdom,Phylum,Class,Order,Family,Genus,Species)%>%   
#  dplyr::summarize(Mean =                      
#                     mean(Abundance, na.rm=TRUE)) %>%    
#  filter(Mean > 0.00) %>%                                
#  arrange(Class)  

## `summarise()` has grouped output by 'Season', 'Kingdom', 'Phylum', 'Class',
## 'Order'. You can override using the `.groups` argument
#Creating a new column that combines both the Class and Family level information #so that Family level identifiers that are not unique to one Class (such as Ambiguous_taxa) aren't merged into a single large category when graphed 
PRSeqR_family_Location$Genus_species <- paste(PRSeqR_family_Location$Class, PRSeqR_family_Location$Species, sep="_")                                                 #Creating plot   
## Warning: Unknown or uninitialised column: `Species`.
PRSeqR_Location_family_bar <- ggplot(PRSeqR_family_Location, aes(x = Location, y = Mean, fill = Genus_species)) +  geom_bar(stat="identity", colour = "black", size=0.3, position = "fill") +   scale_fill_manual(values = c("black",                                                                 "#440154FF", "#450659FF","#460B5EFF",                                                                "#472D7AFF",                                                                 "#3B518BFF", "#3A548CFF", "#38598CFF", "#365C8DFF",                                 "#34608DFF", "#33638DFF", "#31678EFF", "#306A8EFF",                                "#2E6E8EFF", "#2D718EFF", "#2B748EFF",                                "#2A778EFF", "#297B8EFF",                                                                "#1F958BFF", "#1F988BFF", "#1E9C89FF", "#1F9F88FF",                                "#1FA287FF", "#21A585FF", "#23A983FF",                                "#25AC82FF", "#29AF7FFF", "#2DB27DFF",                                "#32B67AFF", "#37B878FF", "#3CBC74FF",                                                                "#57C766FF", "#5EC962FF",                                                                "#A2DA37FF",                                                                "#DAE319FF", "#E4E419FF",  "#ECE51BFF", "#F5E61FFF",                                                                "darkorchid1", "darkorchid2", "darkorchid3",                                                                 "#A21C9AFF", "#A62098FF", "#AB2394FF", "#AE2892FF", "#B22B8FFF", "#B6308BFF", "#BA3388FF", "#BE3885FF", "#C13B82FF", "#C53F7EFF",                                "#C8437BFF", "#CC4678FF", "#CE4B75FF", "#D14E72FF", "#D5536FFF", "#D7566CFF", "#DA5B69FF", "#DD5E66FF", "#E06363FF","#00CCCC", "#006600",                                                                "#FF6699",                                                                "#F68D45FF", "#FCA537FF",                                                                "#F6E726FF", "#F4ED27FF",                                                                                                 "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018",                                                                "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018",                                "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018",                                "#00489C", "#CCCCCC", "#999999", "#A1C299", "#300018","#191970","#353935","#F0FFFF","#89CFF0","#0000FF","#7393B3","#088F8F","#0096FF","#5F9EA0","#0047AB","#6495ED","#00FFFF","#00008B","#6F8FAF","#1434A4","#7DF9FF","#6082B6","#00A36C","#3F00FF","#5D3FD3","#ADD8E6","#191970","#000080","#1F51FF","#A7C7E7","#CCCCFF","#B6D0E2","#96DED1","#4169E1","#0F52BA","#9FE2BF","#87CEEB","#4682B4","#008080","#40E0D0","#0437F2","#40B5AD","#0818A8","#AAFF00","#5F9EA0","#097969","#AFE1AF","#DFFF00","#E4D00A","#00FFFF","#023020","#50C878","#5F8575","#4F7942","#228B22","#7CFC00","#F2D2BD","#FFAC1C","#CD7F32","#DAA06D","#CC5500","#E97451","#E3963E","#F28C28","#D27D2D","#B87333","#FF7F50","#F88379","#DA70D6","#C3B1E1","#CCCCFF","#673147","#A95C68","#800080","#51414F","#953553","#D8BFD8","#630330","#7F00FF","#722F37","#BDB5D5","#FCF55F","#FFFFF0","#F8DE7E","#F0E68C","#FAFA33","#FBEC5D","#F4BB44","#FFDB58")) +      guides(fill = guide_legend(reverse = FALSE, keywidth = 1, keyheight = 1)) +   ylab("Relative Abundance (Species > 0) \n") + xlab("Location") +   theme_minimal()+theme(panel.border = element_rect(fill=NA, colour = "black"),    panel.grid.minor = element_blank(),   panel.grid.major = element_blank(),    axis.text.x  = element_text(angle = 90, vjust = 0.5, hjust=1, size=8,                                colour="black"),    axis.text.y = element_text(size=8, colour="black"),    plot.title = element_text(hjust = 0.5),    axis.ticks.x = element_line(colour="#000000", size=0.1),    axis.ticks.y = element_line(colour="#000000", size=0.1),   legend.text = element_text(size = 7))

print(PRSeqR_Location_family_bar)

p8 <- ggplotly(PRSeqR_Location_family_bar)
#(2/2) This step is needed to create an interactive plot on gtihub. After it has been pushed to github you can download the html doc and it will open as the interactive plot via your web browser.
htmlwidgets::saveWidget(p, "PRSeqR_Location_family_bar")

htmltools::tags$iframe(
  src=file.path(getwd(), "PRSeqR_Location_family_bar"),
  width="100%",
  height="600",
  scrolling="no",
  seamless="seamless",
  frameBorder="0"
)

Transforming Phyloseq object to relative abundance

#Transform to relative abundance  
PR_RA_Normal = transform_sample_counts(PR_Normal, function(x){x / sum(x)})    

#View the otu_table to see the transformed relative abundances  
View(otu_table(PR_RA_Normal)) 

NMDS Plots

Choose colors and shapes

#Choose colors
plot.colors <- c("steelblue2","purple4","darkorange","firebrick","springgreen4", "gold", "darkblue", "yellowgreen","turquoise4", "orange","indianred","darkslategrey", "lightblue","darkgreen","mediumaquamarine","gray48","mediumorchid1", "#5F7FC7","#DA5724", "#508578", "#CBD588","#CD9BCD","#AD6F3B", "#673770","#D14285", "#652926", "#C84248", "#8569D5", "#5E738F","#D1A33D","#8A7C64", "#599861","dodgerblue","darkmagenta", "forestgreen","steelblue1", "cyan","mediumorchid3", "cadetblue3", "yellow", "#00FF7F","#40E0D0","#E97451","#0000FF","#7393B3","#088F8F","#0096FF","#5F9EA0","#0047AB","#6495ED","#00FFFF","#00008B","#6F8FAF","#1434A4","#7DF9FF","#6082B6","#00A36C","#3F00FF","#5D3FD3","#ADD8E6","#191970","#000080","#1F51FF","#A7C7E7","#CCCCFF","#B6D0E2","#96DED1","#4169E1","#0F52BA","#9FE2BF","#87CEEB","#4682B4","#008080","#40E0D0","#0437F2","#40B5AD","#0818A8","#AAFF00","#5F9EA0","#097969","#AFE1AF","#DFFF00","#E4D00A","#00FFFF","#023020","#50C878","#5F8575","#4F7942","#228B22","#7CFC00","#F2D2BD","#FFAC1C","#CD7F32","#DAA06D","#CC5500","#E97451","#E3963E","#F28C28","#D27D2D","#B87333","#FF7F50","#F88379","#DA70D6","#C3B1E1","#CCCCFF","#673147","#A95C68","#800080","#51414F","#953553","#D8BFD8","#630330","#7F00FF","#722F37","#BDB5D5","#FCF55F","#FFFFF0","#F8DE7E","#F0E68C","#FAFA33","#FBEC5D","#F4BB44","#FFDB58", "#C4B454")
plot.shape <- c(21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21)

Subset the samples

Subsetting from PR_physeq_RA_Normal, which includes only the non-enriched samples. This dataframe has been filtered so that ASVs that were 0 in all samples are excluded.

#PR_S003 <- PR_physeq_RA_Normal %>% subset_samples(Project %in% c("S003"))

#PR_S015 <- PR_physeq_RA_Normal %>% subset_samples(Project %in% c("S015"))

#PR_2021 <- PR_Normal %>% subset_samples(Year %in% c("2021"))

#PR_2022 <- PR_Normal %>% subset_samples(Year %in% c("2022"))

NMDS of both runs

all.nmds.source.ord <- ordinate( 
  physeq = PR_RA_Normal, 
  method = "NMDS", 
  distance = "bray")
## Run 0 stress 0.1850681 
## Run 1 stress 0.1850731 
## ... Procrustes: rmse 0.0007928043  max resid 0.00578824 
## ... Similar to previous best
## Run 2 stress 0.2139555 
## Run 3 stress 0.2207744 
## Run 4 stress 0.2022847 
## Run 5 stress 0.2123585 
## Run 6 stress 0.1850805 
## ... Procrustes: rmse 0.001712935  max resid 0.01392611 
## Run 7 stress 0.1871644 
## Run 8 stress 0.1871644 
## Run 9 stress 0.2024291 
## Run 10 stress 0.2035301 
## Run 11 stress 0.2080844 
## Run 12 stress 0.1850833 
## ... Procrustes: rmse 0.003884105  max resid 0.03295702 
## Run 13 stress 0.1850062 
## ... New best solution
## ... Procrustes: rmse 0.005042124  max resid 0.03300655 
## Run 14 stress 0.2179661 
## Run 15 stress 0.2297205 
## Run 16 stress 0.1849911 
## ... New best solution
## ... Procrustes: rmse 0.002372687  max resid 0.01513598 
## Run 17 stress 0.1849877 
## ... New best solution
## ... Procrustes: rmse 0.001707955  max resid 0.015439 
## Run 18 stress 0.1850491 
## ... Procrustes: rmse 0.005735748  max resid 0.05125967 
## Run 19 stress 0.2237442 
## Run 20 stress 0.1870003 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##      1: no. of iterations >= maxit
##     18: stress ratio > sratmax
##      1: scale factor of the gradient < sfgrmin

NMDs plot of location for combined runs

#Plot, color coding by particle type
all.nmds.Location <- plot_ordination(
  physeq = PR_RA_Normal,
  ordination = all.nmds.source.ord) + 
  scale_fill_manual(values = plot.colors, "Location")+
  scale_shape_manual(values = plot.shape, name = "Location") +
  geom_point(mapping = aes(fill = factor(Location), shape = Location, size = 5)) +
  guides(size=FALSE) +
  guides(shape = guide_legend(override.aes = list(size = 3))) +
  theme(plot.title = element_text(size = 18),
        text = element_text(size = 18), 
        axis.title = element_text(size = 15),
        panel.spacing = unit(1, "lines"), 
        panel.border = element_rect(colour = "black", fill = NA, size = 0.5), 
        panel.background = element_blank(), 
        legend.text = element_text(size = 15),
        legend.title = element_text(size = 15),
        legend.position = c(.99, .99),
        legend.justification = c("right", "top"),
        legend.box.just = "right",
        legend.margin = margin(6, 6, 6, 6)) +
  stat_ellipse(aes(group=Location))
## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
print(all.nmds.Location)

NMDS plot of season

#Plot, color coding by seaosn
all.nmds.Season <- plot_ordination(
  physeq = PR_RA_Normal,
  ordination = all.nmds.source.ord) + 
  scale_fill_manual(values = plot.colors, "Season") +
  scale_shape_manual(values = c("D" = 21, "W" = 21), name = "Season") +
  geom_point(mapping = aes(fill = factor(Season), shape = Season, size = 5)) +
  guides(size=FALSE) +
  guides(shape = guide_legend(override.aes = list(size = 3))) +
  theme(plot.title = element_text(size = 18),
        text = element_text(size = 18), 
        axis.title = element_text(size = 15),
        panel.spacing = unit(1, "lines"), 
        panel.border = element_rect(colour = "black", fill = NA, size = 0.5), 
        panel.background = element_blank(), 
        legend.text = element_text(size = 15),
        legend.title = element_text(size = 15),
        legend.position = c(.99, .99),
        legend.justification = c("right", "top"),
        legend.box.just = "right",
        legend.margin = margin(6, 6, 6, 6)) +
  stat_ellipse(aes(group=Season))

print(all.nmds.Season)

PERMANOVA Analysis

#filt_complete_bray <- phyloseq::distance(PR_RA_Normal, method = "bray")
#filt_complete_bray_df <- data.frame(sample_data(PR_RA_Normal))
#filt_complete <- adonis2(filt_complete_bray ~ mm_Hg*DO_mg_L*DO_LDO*Specific_Conductance_uS_cm, data = filt_complete_bray_df)
#filt_complete
#Run PERMANOVA with Bray-Curtis dissimilarity
#filt_complete_bray <- phyloseq::distance(PR_RA_Normal, method = "bray")
#filt_complete_bray_df <- data.frame(sample_data(PR_RA_Normal))
#filt_complete <- adonis2(filt_complete_bray ~ Temperature_C*Location*pH*Salinity_ppt, data = filt_complete_bray_df)
#filt_complete
#Run PERMANOVA with Bray-Curtis dissimilarity
#filt_complete_bray <- phyloseq::distance(PR_RA_Normal, method = "bray")
#filt_complete_bray_df <- data.frame(sample_data(PR_RA_Normal))
#filt_complete <- adonis2(filt_complete_bray ~ Temperature_C*Location*DO_LDO*pH*Salinity_ppt, data = filt_complete_bray_df)
#filt_complete
#Convert output to dataframe
#filt.complete.pvalue <- as.data.frame(filt_complete)

#Create table of p-values
#knitr::kable(filt.complete.pvalue, row.names = NA) %>% 
 # kableExtra::kable_styling("striped", 
                           # latex_options="scale_down") %>% 
  #kableExtra::scroll_box(width = "100%")

PERMANOVA analysis of location vs pH

#Subset phyloseq object
#filt.mp.particle.RA <- subset_samples(MPfiltRA, sample_type == "Particle")
#sample_data(filt.mp.particle.RA)
#Run PERMANOVA with Bray-Curtis dissimilarity
#filt_particle_bray <- phyloseq::distance(filt.mp.particle.RA, method = "bray")
#filt_particle_bray_df <- data.frame(sample_data(filt.mp.particle.RA))
#filt_all <- adonis2(filt_particle_bray ~ particle_type*effluent*week*polymer_type, data = filt_particle_bray_df)
#filt_all
#Convert output to dataframe
#filt.all.pvalue <- as.data.frame(filt_all)

#Create table of p-values
#knitr::kable(filt.all.pvalue, row.names = NA) %>% 
  #kableExtra::kable_styling("striped", 
                          # latex_options="scale_down") %>% 
  #kableExtra::scroll_box(width = "100%")
#Run PERMANOVA with Bray-Curtis dissimilarity, only looking at particle type (glass vs. plastic)
#filt_particle_only_bray <- phyloseq::distance(filt.mp.particle.RA, method = "bray")
#filt_particle_only_bray_df <- data.frame(sample_data(filt.mp.particle.RA))
#filt_only_particle <- adonis2(filt_particle_only_bray ~ particle_type, data = filt_particle_only_bray_df)
#filt_only_particle

NMDS Plot extras in R: Envfit

##Load libraries and read in your data
library(vegan)
library(ggplot2)
##Subset your data to only environmental data and abundance data. 
##I already had a dataframe from the PR_RA_Normal for abundance data, so I just needed to do the same for the sample data (aka environmental data).
#df_env <- as.data.frame(sample_data(PR_RA_Normal))
#df_com <- as.data.frame(otu_table(PR_RA_Normal))
## Individual data frames created, now I am selecting the columns that will be used for this plot and renaming the data frame to com and env. 
#com = df_com[,1:95]
#env = df_env[,15:17, 20, 21,24]
#convert com to a matrix and perform the NMDs ordination
#m_df_com = as.matrix(com)

#nmds code
#set.seed(123)
#nmds = metaMDS(m_df_com, distance = "bray")
#nmds
##Running envfit function with env dataframe
#en = envfit(nmds,env, permutations = 999, na.rm = TRUE)

RDA - Redundancy Analysis

library(phyloseq)
library(ggplot2)
library(microViz)

# install.packages("remotes")
#remotes::install_github("statdivlab/corncob")
#library(corncob)
#> microViz version 0.12.0 - Copyright (C) 2023 David Barnett
#> ! Website: https://david-barnett.github.io/microViz
#> ✔ Useful?  For citation details, run: `citation("microViz")`
#> ✖ Silence? `suppressPackageStartupMessages(library(microViz))`
#knitr::opts_chunk$set(fig.width = 7, fig.height = 6)
#PR_RA_Normal <- corncob::PR_phylo
#PR_phylo
#> phyloseq-class experiment-level object
#> otu_table()   OTU Table:         [ 36349 taxa and 91 samples ]
#> sample_data() Sample Data:       [ 91 samples by 15 sample variables ]
#> tax_table()   Taxonomy Table:    [ 36349 taxa by 7 taxonomic ranks ]
##Redundancy analysis is a constrained ordination method. It displays the microbial variation that can also be explained by selected constraint variables. Behind the scenes, a linear regression model is created for each microbial abundance variable (using the constraints as the explanatory variables) and a PCA is performed using the fitted values of the microbial abundances.

#PR_RA_Normal %>%
 # ps_mutate(
  #  IBD = as.numeric(ibd == "ibd"),
   # Female = as.numeric(gender == "female"),
   #  ) %>%
 # tax_transform("clr", rank = "Genus") %>%
 # ord_calc(
    #constraints = c("IBD", "Female", "Abx."),
    # method = "RDA", # Note: you can specify RDA explicitly, and it is good practice to do so, but microViz can guess automatically that you want an RDA here (helpful if you don't remember the name?)
  #  scale_cc = FALSE # doesn't make a difference
 # ) %>%
 # ord_plot(
 #   colour = "DiseaseState", size = 2, alpha = 0.5, shape = "active",
 #   plot_taxa = 1:8
 # )